#include "aggAux.h"
#include "manager.h"

// Performs the tests of the control volumes and edges coefficients, to
// check if the agglomeration procedure was successfully
// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
// IN: aggMesh[iMesh], setOfDomains
// OUT: Number of CVs with errors

void aggTest (pMesh theMesh, set<int> &setOfDomains)
{
// Definig the domain flag
set<int>::iterator iter = setOfDomains.begin();

// Constant defining tolerance
const double TOL(1.0E-10);
// Variable to store the number of edges with false coeff.
int nErrorCV(0);
// Variables used to store the coefficients
double val, deltaDij;
char *Cij_string[2], *Dij_string[2], *deltaDij_string[2];
for (int i=0; i<2; i++)
{
	Cij_string[i] = new char[64];
	Dij_string[i] = new char[64];
	deltaDij_string[i] = new char[64];
}

vector<pVertex> vertex;

// Coefficients test (Sum(C's) + Sum(D's) for each node = 0)
cout<< "   Testing Coefficients...";
VIter vit = M_vertexIter(theMesh);
while (pEntity node = VIter_next(vit)){

	int nodeID = EN_id(node);
	double sum[2] = {.0,.0};

	int numEdges = V_numEdges(node);        // get number of edges surrounding
                                            // the entity vertex
    for (int i=0; i<numEdges; i++)
    {

        pEdge edge = V_edge(node,i);
        if ( !theMesh->getRefinementDepth(edge) ){

        M_GetVertices(edge,vertex);
        int id0 = EN_id(vertex[0]), id1 = EN_id(vertex[1]);
		double signal = (nodeID == id1)?-1.0:1.0;
		std::vector<double> Cij(3,.0),Dij(3,.0);

        // Summing C's coefficients (Internal)
        sprintf(Cij_string[0],"Cx-%d",*iter);
        sprintf(Cij_string[1],"Cy-%d",*iter);
        for (int j=0; j<2; j++)
        {
            val = .0;
            EN_getDataDbl(edge,MD_lookupMeshDataId(Cij_string[j]),&val);
            sum[j]+= signal*val;
        }

        // Summing D's coefficients (Boundaries)
        sprintf(Dij_string[0],"Dx-%d",*iter);
        sprintf(Dij_string[1],"Dy-%d",*iter);
        sprintf(deltaDij_string[0],"deltaDx-%d",*iter);
        sprintf(deltaDij_string[1],"deltaDy-%d",*iter);
        for (int j=0; j<2; j++)
        {
            val = .0;
            deltaDij = .0;
            EN_getDataDbl(edge, MD_lookupMeshDataId(Dij_string[j]),&val);
            EN_getDataDbl(edge, MD_lookupMeshDataId(deltaDij_string[j]),&deltaDij);
            //if (val !=.0) cout<<"D ["<<nodeID<<"]: "<<val<<endl;
            if (nodeID == id0)
            {
                sum[j] = sum[j] + val + deltaDij;
            }else
            {
                sum[j] = sum[j] + val - deltaDij;
            }
        }
        //printf("[%d %d] Cij: %E %E\tsignal: %f\tDij: %E %E\n",id0,id1,Cij[0],Cij[1],signal,Dij[0],Dij[1]);
        vertex.clear();
        }
    }
    if (fabs(sum[0]) > TOL || fabs(sum[1]) > TOL)
    {
        printf("   [%d] - Cij+Dij = %E %E\n",nodeID,sum[0],sum[1]);//throw 1;
        nErrorCV++;
    }
}
VIter_delete(vit);//throw 1;
// Printing the total number of Control Volumes with false coefficients (sum(C) + sum(D) != 0)
cout << "   Number of CVs with error: " << nErrorCV << "\n";

// Total volume test(Sum(Control Volumes) = total volume)
cout<< "   Testing Volumes...";
double s = .0;
vit = M_vertexIter(theMesh);
while (pEntity node = VIter_next(vit)){
	int nodeID = EN_id(node);
	char volume_string[64];
	sprintf(volume_string,"volume-%d",*iter);
   	double volume = .0;
   	EN_getDataDbl(node,MD_lookupMeshDataId(volume_string),&volume);
	s += volume;
}
VIter_delete(vit);//throw 1;
printf("        Volume total = %E\n",s);
}
